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ABSTRACT 

This series of papers investigates the early stages of planet formation by modeling the evolution of the gas 
and solid content of protostellar disks from the early T Tauri phase until complete dispersal of the gas. In this 
first paper, I present a new set of simplified equations modeling the growth and migration of various species of 
grains in a gaseous protostellar disk evolving as a result of the combined effects of viscous accretion and photo- 
evaporation from the central star. Using the assumption that the grain size distribution function always maintains 
a power-law structure approximating the average outcome of the exact coagulation/shattering equation, the model 
focuses on the calculation of the growth rate of the largest grains only. The coupled evolution equations for 
the maximum grain size, the surface density of the gas and the surface density of solids are then presented and 
solved self-consistently using a standard l-i-l dimensional formalism. I show that the global evolution of solids 
is controlled by a leaky reservoir of small grains at large radii, and propose an empirically derived evolution 
equation for the total mass of solids, which can be used to estimate the total heavy element retention efficiency 
in the planet formation paradigm. Consistency with observation of the total mass of solids in the Minimum Solar 
Nebula augmented with the mass of the Oort cloud sets strong upper limit on the initial grain size distribution, as 
well as on the turbulent parameter at. Detailed comparisons with SED observations are presented in a following 
paper 

Subject headings: accretion disks - methods: numerical - solar system: formation 
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I. INTRODUCTION 
1.1. Theoretical and observational motivations 

More than two hundred extrasolar planets have now been 
detected, revealing surprising diversity. Doppler surveys have 
provided a large database of masses, orbital radii and eccentric- 
ities, which show notably few (and a few notable) systematics, 
as for example the relationship between stellar metallicity and 
the number of detected planets (Fischer & Valenti, 2005). Tran- 
sit detections are now also beginning to show a large diversity 
in the internal structure of planets with otherwise very similar 
properties (Guillot et ai, 2006). 

Fast-forwarding back a few Gyr, one can rightfully expect to 
find the origin of exo-planetary diversity in the equivalent diver- 
sity of protostellar disks. And evidence has indeed been found 
to support this idea. The observed fraction of stars showing 
excess at near-IR (Haisch et al. 2001, Hartmann et al. 2005, 
Sicilia-Aguilar et al. 2006) and/or mid-IR wavelengths (Ma- 
majek et al. 2004) steadily decreases from nearly 100% for 
stars within the youngest clusters, to zero for stars within clus- 
ters older than about 20 Myr. This correlation has long been 
interpreted as clear evidence for disk dispersal within a typical 
timescale of about lOMyr, but is now beginning to gather addi- 
tional interest as evidence for a large variation in the disk dis- 
persal rates amongst similar type stars within the same cluster. 
This dispersion could be related to variations in the initial disk 
conditions and/or to the characteristics of the host star (Hueso & 
Guillot, 2005). Other possible tracers of disk structure and/or 
evolution (such as the crystallinity fraction and grain growth) 
also reveal significant diversity: for instance, co-eval stars of 
similar types show evidence for very different crystallinity frac- 
tions (Meeus et al. 2003 for T Tauri stars, Apai et al. 2005 for 
brown dwarves). 



Can the origin of this dynamical and structural diversity in- 
deed be traced back to the initial conditions of the disk? Quali- 
tatively speaking, can it explain why some systems form planets 
while others don't? Quantitatively speaking, is there a link be- 
tween the initial angular momentum and mass of the disk and 
the characteristics of the emerging planetary system? 

Meanwhile, stringent upper bounds on the total amount 
of heavy elements typically remaining as planetary building 
blocks have been deduced from the very low metallicity dis- 
persion measured amongst similar type stars within the same 
cluster by Wilden et al. (2002). This result is puzzling in the 
light of the contrastingly large range of observed disk survival 
timescales: how can widely different dynamics lead to similar 
retention efficiencies. 

A necessary step towards answering these questions is the de- 
velopment of a comprehensive numerical model capable of fol- 
lowing the formation and evolution of planetary systems from 
their earliest stages to the present day, including all of the phys- 
ical processes currently understood to play a role in the evolu- 
tion of the gas and solids. 

The standard core-accretion model of planet formation be- 
gins with the condensation of heavy elements into small grains, 
followed by their stochastic collisional growth into successively 
larger aggregates until they reach a typical mass (either col- 
lectively or individually) where mutually induced gravitational 
forces begin to influence their motions. The small planetes- 
imals then continue growing by accreting each other (together 
with some of the disk gas), until a critical point is reached where 
runaway gas accretion may eventually begin. This first plane- 
tary formation phase ends with the dispersal of the disk gas, 
possibly by photo-evaporation, although gravitational interac- 
tions between the various bodies continue taking place result- 
ing in close encounters (sometimes collisions) with dynamical 
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rearrangement of the system (including ejection, shattering, co- 
agulation). 

In this paper I present a numerical model for the first stage of 
this process, in which a protostellar disk and all of its contents 
(both in gaseous and in solid form) are evolved simultaneously 
until complete dispersal of the gas. The next stages of evolu- 
tion from this point onward are best treated with an N-body 
code, for which the results presented here could be used as ini- 
tial conditions. 

Recent data obtained with the Spitzer Space Telescope has 
provided valuable information on the evolution of grains in pro- 
tostellar disks, which can be used to both construct and test the 
desired planet formation model. Since the near- and mid-IR 
ranges of the observed spectral energy distributions (SEDs) are 
essentially due to reprocessing of the stellar radiation by small 
dust grains, the key to modeling planet formation in the con- 
text of evolving disks is to better understand the relationship 
between the observable SEDs and the physics which couple the 
gas and dust dynamics under the gravitational and radiative in- 
fluence of a central star. This is done in Paper II (Alexander & 
Garaud, 2007). 

1.2. General methodology 

This work presents a new versatile numerical tool to study 
the evolution of both gas and solids in protostellar disks, from 
classical T Tauri disks to transition disks and finally to forming 
planetary systems (embedded perhaps in a debris disk). The 
model developed takes into account the following physical phe- 
nomena: (i) axisymmetric 1h-1D gas dynamics around the cen- 
tral star, (ii) photo-evaporation by the central star, (iii) continu- 
ous grain size distribution maintained by growth and fragmen- 
tation, (iv) grain sublimation and condensation, (v) multiple 
grain species (iron, silicates, ices), (vi) gas-grain coupling in- 
cluding turbulent dust suspension, turbulent diffusion and drift 
and (vii) gravitational interaction between forming embryos (in 
a statistical sense). 

While the general goal of modeling the early disk evolution 
has been pursued by many others before, this particular model 
is the first to include all of the physics listed above in a single, 
well-tested, fast and practical algorithm. Other physical phe- 
nomena such as photo-evaporation by nearby stars, truncation 
of the disk by stellar fly-by, or planetary migration are easy to 
implement, but not discussed here. In order to place the model 
in context, it is useful to summarize briefly existing work on the 
subject. A more thorough discussion of the results in the light 
of previous work can be found in ^ 

Axisymmetric gas dynamics in a viscously dominated ac- 
cretion disk has been thoroughly analyzed by Lynden-Bell & 
Pringle (1974). In subsequent work, particular attention was 
given to studying the disk structure and evolution in the light 
of SED observations (see Hartmann et al. 1998 for example). 
Photo-evaporation of the gas by UV photons (either ambient 
and/or emerging from central star) is now thought to play a 
major role in the dispersal of the disk gas. This was studied 
in detail by Hollenbach et al. (1994), and later proposed by 
Clarke, Gendrin & Sotomayor (2001) as a possible model pro- 
viding the characteristic "two-timescale" evolution (namely a 
long lifetime with a rapid dispersal time) required by the low 
relative abundance of transition disks (see the reviews by Hol- 
lenbach & Gorti 2005, and Dullemond et al. 2007). 

Meanwhile, the study of the evolution of solids in protostel- 
lar disks also has a long history, where the particular emphasis 



has in the vast majority of cases been to model the formation 
of our own solar system. The early works of Whipple (1972) 
and Weidenschilling (1977) laid the foundation for studying the 
motion of small solid bodies in the early solar nebula. VoeUc et 
al. (1980) developed a theory for the dynamical coupling of 
solid particles with turbulent eddies, which enabled many fur- 
ther studies of the collisional growth of dust grains into plan- 
etesimals (Weidenschilling, 1984 and subsequent papers, Wei- 
denschilling & Cuzzi 1993, Stepinski & Valageas 1997, Suttner 
& Yorke 2001, Dullemond & Dominik, 2005). Finally, steady 
progress in the interpretation of various cosmochemistry data 
has prompted the need for a better understanding of the evo- 
lution of the various chemical species present in the disk, in 
particular water In addition to their own work, Ciesla & Cuzzi 
(2006) present an excellent review of recent advances in the 
field. 

Combining the evolution of solids with the evolution of the 
gas with the aim of bridging the gap between SED interpreta- 
tions and our own solar system formation is naturally the next 
step in this scientific exploration process. The work of Sut- 
tner & Yorke (2001) pioneered the concept when looking at 
grain growth and migration in the very early stages of the disk 
formation (first few 10"^ yr). Alexander & Armitage (2007) 
(AA07 hereafter) were recently the first to combine state-of- 
the-art photo-evaporation models with grain migration to gain 
a better understanding of the nature of some forming transi- 
tion disks. The proposed model draws from many of the fun- 
damental ideas of these previous studies; in particular, it can 
be thought of as a generalization of the AA07 model which 
includes the effects of grain growth, sublimation and condensa- 
tion. 

Theoretical studies of dust growth typically require the so- 
lution of a collisional equation at every spatial position of the 
disk. Amongst some of the difficulties encountered one could 
mention the determination of the particle structure, the sticking 
efficiency, the shattering threshold and the size distribution of 
the fragments, and not least the relative velocities of the parti- 
cles before collision. Indeed, while the motion of particles in a 
laminar disk is fairly easy to compute, matters are complicated 
when dynamical coupling between grains and turbulent eddies 
is taken into account. Tiny grains are well-coupled with the 
gas though frictional drag, while larger "boulders" only feel the 
eddies as a random stochastic forcing. The intrinsic dispersion 
and the relative velocities of the particles can be modeled statis- 
tically provided one assumes the gas eddies follow a turbulent 
Kolmogorov cascade from the macro-scale to the dissipation 
scale. This idea was originally proposed by Voelk et al. (1980) 
and more recently reviewed by various authors, notably Wei- 
denschiUing (1984). Yorke & Suttner (2001) and DuUemond 
& Dominik (2005) used these velocity prescriptions to evaluate 
the rate of growth of particles in protostellar disks by solving 
the full coagulation equation. Their results show that the colli- 
sional growth of particles in the inner regions of the disk is too 
fast, unless shattering is taken into account. It is therefore vital 
to include it in evolutionary models of disks as well. 

However, solving for the complete coagulation/shattering 
equation for every particle size, at every timestep and for every 
position in the disk is computationally prohibitive. Statistical 
surveys of the typical outcome of the disk evolution for a wide 
range of stellar parameters and initial conditions cannot be done 
in this fashion. 

The novel part of this work concerns the modeling of the evo- 
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lution of the grain size distribution function under collisional 
coagulation and shattering. The underlying assumption of the 
model proposed is that collisions between dust grains are fre- 
quent enough for a quasi-steady coagulation/shattering balance 
to be achieved in such a way as to maintain a power-law parti- 
cle size distribution function with index -3.5 as in the ISM, but 
with varying upper size cutoff .Smax- With this assumption, the 
study of the evolution of solids in the disk can be reduced to 
a small set of one-dimensional partial differential equations for 
the maximum particle size SmaxCr t), the total surface density of 
gas E(r,f), as well as the total surface density of solids and va- 
por for each species considered (I]p(r,f) and Sy(r, f), where / is 
the index referencing the species). Here r is the radial distance 
from the central star and t is time. This idea is to be considered 
as an alternative approach to the work of Ciesla & Cuzzi (2006) 
for instance, who equivalently model the evolution of gas and 
solids in the disk over the course of several Myr, simplifying 
the collision/shattering balance by considering only four "size" 
bins (vapor, grains, rapidly drifting "migrators" and finally very 
large planetesimals). 

1.3. Outline of the paper 

The derivation of the model is presented in complete detail 
in ^(the result-minded reader may prefer to jump straight to 
^and ® . The standard gas dynamics equations together with 
the photo-evaporation model used are well-known, and summa- 
rized for completeness in ^2.1| and 32.21 The basic assumptions 
for the particle size distribution model considered as the ba- 
sis for this paper are presented in 32.31 The stochastic motion 
of solids in the nebula resulting from frictional coupling with 
turbulent eddies and from mutual gravitational encounter have 
been studied by many others before. Key results from these 
works are presented in 32.41 and later used in 32.51 and 32.61 
to derive new equations for the growth of grains into planetes- 
imals, as well as the evolution of the total surface density of 
particles. Finally, 32.71 summarizes the very simple sublima- 
tion/condensation model used here. 

A general overview of the typical inputs and outputs of the 
numerical model are given in |3]and ^respectively. In order to 
gain a better understanding of the numerical results, §5 presents 
existing and new analytical work characterizing the global fea- 
tures of the model (gas dynamics in 35.11 grain growth in 35.21 
evolution of solids in 35.31 §5.4 and §5.5). In particular, a plau- 
sible new semi-analytical evolution equation for the total mass 
of solids in the disk is presented in 35.3.21 which depends only 
on the initial conditions of the disk. Finally, the model and re- 
sults are discussed in ^ Although this paper focuses primarily 
on presenting the methods used (while paper II discusses the 
observable properties of the modeled disks), I give some esti- 
mates for the heavy-element retention efficiency of disks as a 
function of the model parameters, and show how one could rec- 
oncile the high diversity of observed disk properties with the 
low dispersion in metallicities for star within the same cluster 
(Wilden et al. 2002). Conclusions are summarized in ^ 

2. MODEL SETUP 

2. 1 . Evolution of the gas disk 

In all that follows, I assume that the gas disk evolves inde- 
pendently of the solids. Note that this is only true as long as the 
surface density of the gas is much larger than the surface den- 
sity of solids; when the metallicity Z(r,t) = Sp/E approaches 
or exceeds unity, solids begin to influence the evolution of the 



gas through angular momentum exchange and possible gravita- 
tional instabilities. Barring these cases, the standard evolution 
equation for I](r, t) is 

^ + i|.(™S) = -Ew, (1) 
ot r or 

where u is the typical radial velocity of the gas required by con- 
servation of angular momentum in the accretion disk, 

and Ew (where the dot from here on always denotes differentia- 
tion with respect to the time t) is the gas photo -evaporation rate 
modeled following the parametrization of AA07 (see Appendix 
A). 

The gas turbulent diffusivity Vi is modeled using the standard 
a -model 



Vi = Uich = Ui^J^^-^h 



(3) 



where c is the local sound speed and 7 is the adiabatic index of 
the gas. Note that there is a degeneracy between models with 
constant a, and one particular temperature profile, and models 
with non-constant at and another temperature profile yielding 
the same value of Vt. This degeneracy combined with the crude 
Qf-parametrization of turbulent transport used justifies the se- 
lection of a very simple temperature profile: 

T„(r)=Trl^, (4) 

where tau is the distance to the central star in astronomical 
units. The scaleheight of the disk then varies as 



hir) = -hr%'^^' 



(5) 



In what follows, I adopt the same disk model as that used by 
AA07: 

^ = -1/2, 

W = 0.0333. (6) 

Note that AA07 define q as the power index of h(r) instead 
of the power index of T„,{r) used here; the apparently differ- 
ent values do correctly represent the same model. Although 
the numerical algorithm I have developed can be used with any 
input for q and /iau, this particular value of q is preferred as 
it greatly simplifies the analytical interpretation of the numeri- 
cal results; indeed, in this case ft scales linearly with radius, a 
feature which turns out to be particularly useful. 

2.2. Evolution of vapor species 

Chemical species in vapor form are evolved separately using 
the following standard advection-diffusion equation for a con- 
taminant in a fluid of density E moving with velocity u: 



9Et I d ^ 1 d 

ot r or r or 



rviY, 



(7) 



where it was implicitly assumed that the diffusivities of each 
chemical species are equal to the gas viscosity, and u is given 
by equation (|2|i. Sublimation and condensation are assumed to 
be instantaneous on the timescales considered and are calcu- 
lated as a separate numerical step (see 32.7l i. 

2.3. Particle size distribution function 

Collisional encounters between solid particles can result in 
their coagulation or mutual shattering, the latter sometimes fol- 
lowed by the re-accretion of material onto the largest remaining 
fragments. However complex the mechanisms considered are. 



4 



the size distribution function of the particles is naturally ex- 
pected to relax to a quasi-steady equilibrium power-law within 
a few colUsion times. Theoretical arguments on the steady-state 
nature of the coUisional cascade imply that the power-law in- 
dex depends on the relationship between the relative velocities 
of the objects and their material strengths (O'Brien & Green- 
berg, 2003). Such power-laws are observed in the ISM (with 
index -3.5, Mathis, Rumpl & Nordsieck, 1977), for Kuiper-belt 
objects (with varying index depending on the size range) and 
for asteroid-belt objects. This model is constructed by assum- 
ing that encounters are frequent enough to maintain a quasi- 
steady equilibrium, which results in a power-law size distribu- 
tion (with fixed index -3.5) for all particles of size less than 
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d^ 



•Smax 



-3.5 



for 5* G [5, 
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dn 

— = otherwise 
as 



(8) 



where I allow the normaUzing density rimax, and the maximum 
particle size 5inax Vary both with radius and with time. The 
minimum particle size Smm is fixed, although its value does 
not influence the dynamical evolution of the disk as long as 
■Smax ^ '^min (siucc most of the solid mass is contained in the 
largest grains). Note the value of Smm influences the SED since 
the smallest grains contribute the most to the total emitting sur- 
face area. 

If the particles are spherical with uniform solid density ps 
then the total density of solids is 



„1,X 

— m(s)ds = ZHmaxOTmax 

as 



(9) 



provided imin ^ -Smax. where m{s) is the mass of particles of size 
s, and mmax is the mass of particles of size Smax namely 



4n 



m„ 



(10) 



This power-law size distribution function implies that 50% 
of the total mass is contained in particles of size s E 



[0.25.?, 



max: ^*max 



The total surface density of particles is 



(11) 



All condensed heavy elements present at a particular radius r 
are assumed to be fully mixed, or in other words, each particle 
has a mixed chemical composition that can vary depending on 
its radial position within the disk. Within this assumption, tiiaax 
can be related to the total density of solids only, and within the 
particle disk (near the disk midplane), is directly related to the 
total surface density of particles via the equation 



(12) 



2WmaxV27r/Zp 



(assuming pp has a Gaussian profile across the disk with scale- 
height hp). Note that the particle scaleheight hp depends on the 
mechanism exciting the intrinsic particle dispersion, which can 
be frictional coupling with turbulent eddies or mutual gravita- 
tional interactions. It is naturally independent of the particle 
species considered. Explicit expressions for hp in these two 
limits are given below. 



2.4. Particle motion 

Motion of particles within the disk can be induced by various 
possible forces: Brownian motion, motion induced by frictional 
drag with the gas and motion induced by interactions with the 
gravitational potential of the central star or that of other large 
planetesimals. The dominant term depends on the particle size. 

Since the only particles considered here have size i^ax, 
Brownian motion is typicaUy negligible. In a turbulent neb- 
ula, particles of various sizes couple via gas drag to the tur- 
bulent eddies and can acquire significant velocities when their 
typical stopping time is comparable with the eddy turnover 
time. Larger particles are only weakly coupled with the gas 
but undergo significant gravitational interactions with each 
other which constantly excite their eccentricities and incUna- 
tions. These mechanisms can be thought of as various kinds 
of stochastic forcing. Finally, non-stochastic forces arise from 
the gravitational potential of the central star, and when com- 
bined with gas drag, can cause particles to sediment towards 
the mid-plane of the disk as weU as spiral inward (occasionaUy 
outward). 

These regimes are now described in more detail. 

2.4. 1 . Turbulence-induced dynamics 

In this section, I summarize existing results on the statistical 
properties of the dust dynamics resulting from their frictional 
coupling with turbulent eddies, and apply them to the problem 
at hand. 

1. Frictional drag. Particles are coupled to the gas through 
frictional drag. The amplitude of the drag force is statisticaUy 
proportional to the relative velocity between the particle and the 
gas, with a proportionality constant that depends on whether the 
particle size is smaller or larger than the mean-free-path of the 
gas molecules Amip (Whipple, 1972). 

If the particle is much smaller than A^fp (Epstein regime), 
drag forces originate from random collisions with the gas 
molecules, and the typical timescale within which the particle 
wiU stop relative to the gas is 



t{s) = . 

pc 



(13) 



If the particle size is much larger than A^fp (Stokes regime) 
then the gas drag is principally caused by the turbulent wake in- 
duced by the particles as it passes through the gas. In this case, 
the particle stopping time is 

SPi 



T{S) = 



(14) 



where a is the typical velocity of the particle with respect to the 
gas, and the constant Co — 0.165 (see Whipple 1972, Garaud, 
Barriere-Fouchet & Lin 2004). 

In what follows, it is useful to define St{s) as the ratio of the 
local stopping time to the local orbital time = 2-k/Q.y^ (Wei- 
denschilhng, 1977), also called the Stokes number: 



Stis)='-^ 

Td 



(15) 



Note that the Stokes number is equally as often defined as 
rJKT(s) by other authors (DuUemond & Dominik 2005 for in- 
stance). 

2. Relative velocities of particles. As first estimated by Voelk et 
al. (1980) and summarized by DuUemond & Dominik (2005) 
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(see also Weidenschilling 1984), particles of various sizes can 
acquire significant relative velocities through their frictional 
coupling with turbulent eddies. This effect depends on the 
relative values of the eddy turnover time and of the particle 
stopping time. For Kolmogorov turbulence with large-scale 
eddy velocity Vp ~ yfoxc and large-scale turnover time compa- 
rable with the dynamical timescale Td, the Reynolds number 
Re = VeT^ / v determines the eddy turnover time at the dissipa- 
tion scale as = Td/?e"'/^. Then, for two particles of respective 
stopping times t{s) and t(s') 



Av(s,s') = 
Av(s,s') = 
Av(s,s') = 

Av(s,s') = 



(T(s)-T(s')f 



1/2 



Td(r(5) + T(i')) 
Ve if t(s') < Td < t(s) 



if t(s'),t(s)<t. 



(16) 



Td + t(s) 

3 

t(s) + t(s') 



1/2 



Td + r(i') 
max(T(s),r(s'))^ 



V'e if Td < t(s'),t(s) . 
1/2 



Td 



Ve Otherwise. 



Note that in the first limit 1 have set y^lnRe/2 = 1 for sim- 
plicity, which is underestimating the true collisional velocity 
by a factor of no more than about 4. This factor will be com- 
pensated for later (see ^2.5\ . Also note that the expression for 
the relative velocities in ( fTTI ) has been corrected from that of 
WeidenschilHng (1984) or Dullemond & Dominik (2005) to 
account for an error pointed out by Ormel & Cuzzi (2007). 

3. Particle diffusion and effective Schmidt number. The stan- 
dard parametrization for the stochastic motion of particles of 
single size s coupled by gas drag to turbulent eddies is through 
the introduction of a turbulent diffusive mass flux f^is) in the 
particle continuity equation, typically 



/t(^)=-pDp(i)V 



P 



(17) 



where the turbulent diffusivity Dp(s) is related to I't through the 
size-dependent Schmidt number 



DM = 



1^1 
Sc{s) 



(18) 



The smallest particles are fully coupled with the gas so that 
Sc(s) ~ 1 if t(s) ^ Td- The standard parametrization for the 
Schmidt number in the case of large particles has long been 
Sc(s) ~ St{s) (see for instance Dubrulle, Morfill & Sterzik, 
1995), so that Sc{s) can be crudely approximated as Sc(s) = 
1 + St{s). Recent numerical and analytical work have shed 
doubts on this formula in favor of Sc(s) cx St^{s) for large par- 
ticles (Carballido, Fromang & Papaloizou 2006) and have also 
questioned the validity of equation ( fTTI l in favor of a different 
formalism involving the equilibrium solution a Fokker-Plank 
equation. Since these very recent studies have not yet been fully 
completed (in particular, they only consider particle diffusion in 
the z-direction and do not propose an alternative formalism for 
the radial diffusion of particles in the disk), 1 continue for the 
moment to adopt the standard parametrization of the Schmidt 
number Sc(s) = l+St(s). 

For a fluid containing a size distribution of particles, the local 
diffusive mass flux of particles is obtained by integrating f^(s) 
across all sizes, yielding 



/t = -OppV(^ 



(19) 



with Dp = i/t / Scsff and the effective Schmidt number Sc^ff being 

(20) 



5Ceff = 



arctan(V5fmax) 
Note that 5ceff is of order unity when Stj 

jc/tt if St, 



0, as expected, 
CO. This is quite different 



while SCf-ii ~ 2^/St^-i 
from the single particle size case, where the Schmidt number 
scales linearly with particle size instead of with y/Smux in the 
decoupled limit. This reflects the fact that smaller particles 
remain well-coupled with the gas even when particles of size 
imax are fully decoupled. 

4. Dust disk scaleheight. Following the work of Dubrulle, 
Morfill & Sterzik (1995), the dust disk scaleheight /ip can be es- 
timated by seeking stationary solutions of the settling/diffusion 
equation 



d 



dz 



(21) 



^-3^(^f^Kr(w)Pp) = ^ 

where the factor of 1/3 arises from the mass-weighted integral 
of the settling velocities over the dust-size distribution function. 
Integrating this equation with height above the disk and assum- 
ing steady-state yields 



hp = h 



1 + - 



2T7 StnTi^xSC, 



eff 



3 atA/7 
where h is the gas scaleheight. 



-1/2 



(22) 



2.4.2. Gravitationally-induced motions 

As described by Kokubo & Ida (2002), the typical velocity 
dispersion of a swarm of planetesimals (which is also equal to 
their typical relative velocities) can be deduced from the bal- 
ance between gravitational excitation by the largest bodies, and 
damping by gas drag. The typical timescale for the excitation 
of the dispersion <j(s) of planetesimals of size s by protoplanets 
of size imax is given by equation (9) of Kokubo & Ida (2002) 

_ 4r^b<i\s)>y^ a(sf 



1^ . 

max 



where In A is the Coulomb logarithm, typically of the order of 
a few (here, I set InA = 3). The typical orbital separation b 
of the emerging protoplanets is of the order of a few Hill radii 
(Kokubo & Ida 2002): 



b = brv,= \<d 



3M, 



1/3 



(24) 



where = 10. The average incUnation of the planetesimals 
< {^{s) >'/^ is assumed to be of the order of the average eccen- 
tricity, so that < e-{s) >'/^= 2 < i-{s) >'/^. Finally, the random 
velocity of the planetesimals is also assumed to be related to 
their average eccentricity by 

G(s)=<e\s)>'l^vy^. (25) 
The timescale for damping of the typical inclination and eccen- 
tricity of the planetesimals is dictated by Stokes drag, namely 

2m(i) 



7dp 



(26) 



^ CuTTs^pais) 

Equating the two timescales yields the velocity dispersion for 
planetesimals of size s in the presence of protoplanets of size 



-(^)=(2 



1/15 



41nA , 27r h 
^St(s)-^- 
3 Cob r 



1/5 



1/3 



(27) 
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As Kokubo & Ida (2002) found, this expression is only weakly 
dependent on the planetesimal size. If the gravitational pertur- 
bations are assumed to be statistically independent, then the rel- 
ative velocities of the planetesimals are equal to their velocity 
dispersion. The weak dependence on size then implies that one 
can approximate the typical scaleheight of the planetesimals as 

h^~<i\s^,,)>'l^r . (28) 

2.5. Particle growth 

In the proposed model, the particle size distribution func- 
tion is parametrized with the power-law form given in equa- 
tion (|8]l, under the assumption that such power-law is natu- 
rally maintained as the quasi-steady state outcome of a coag- 
ulation/shattering balance. The normalization factor Wmax is di- 
rectly related to the total surface density of the dust Sp, while 
the maximum achievable size ^max slowly grows in time as a 
result of occasionally successful coagulation events. 

Following this idea, I model the evolution equation for ^max 
from the standard coagulation equation 



dm„ 



df 



dn 



(i')m(/) Ay( w , ^')A( w , s')eAs' (29) 



where Av(5tnax,'5') is the average relative velocity between par- 
ticles of size 5,nax ^nd size s' , A(.s,„ax, -^0 is the collisional cross- 
section of the two particles and e is the sticking probability 
of the two particles after the collision, or can be alternatively 
thought of as the average mass fraction of the impactor that 
sticks to the target after each collision. Note that in principle e 
could depend on the collisional velocity, on the structure of the 
particles and on their size. In what follows, the function e will 
be chosen to be constant across all sizes and relative velocities 
for simplicity. This approximation is rather unsatisfactory, but 
merely mirrors insufficient knowledge about the exact charac- 
teristics of the dust or larger particles. It can also be thought of 
as a weighted average of the true collisional efficiency across 
all size ranges and all possible impact velocities. 

2.5.1. Growth of particles in the turbulent regime 

For solid particles typically smaller than a few kilometers 
gravitational focusing is negligible (see below). Within this ap- 
proximation, the collisional cross-section of two particles is re- 
duced to the combined geometrical cross-section: 



A{S,S') = TT(S + S'f 



(30) 



Using the expressions derived in §2.4. 1 for the relative ve- 
locities and the particle disk scaleheight, it is now possible to 
re-write equation ( |29] l in a much simpler form. Three limits 
must first be considered: r(5niax) <C Tj,, < r(5niax) < and 

Td < T(i,nax)- 

Case 1: T(imax) <C t^. In this case the particle growth is gov- 
erned by 



dSmax Sp /- h I — 1\ 

df 8/9s /ip 



where the integral l\ is given by 
»i 



'1 



3/2/ 



(31) 



(32) 



Assuming that the sticking efficiency e is constant, and that 
■Smin/^max ^ 1 the integral simpUfies to I\ ~ l.Se. 



Case 2: < r(imax) < ^d. In this case, 

d^max 



df 



where the integral h is given by 
h 



(33) 



(34) 



Under the same assumptions as in Case l,h — 8e. 

Case 3: Td <C T(imax)- This third case is slightly more com- 
plex, as the integral over particle sizes must be split between 
two bins, namely t(s') < Td and t(s') > t^. This yields (in the 
limit considered) 



ds„ 



df 




(35) 



where /3 ~ 2e and I4 ~ 5e5f~ax'- 

For simplicity, the three cases can be combined into one for- 
mula only, namely 



d-^max _ 

df p. 



atSt„ 



'/2pV l+645f2ax(2 + 55f-0ax')-2rd 



(36) 



This expression overestimates the growth rate of the smallest 
particles (i.e. case 1) by a factor of about four This error 
closely compensate for the factor of 4 underestimate in the 
collisional velocity of the smallest particles deliberately made 
in equation ST% . The proposed expression recovers the for- 
mula for grain growth proposed by Stepinski & Valageas (1997) 
within factors of order unity (see their equation (38)). 

2.5.2. Growth of particles in the gravitationally dominated 

regime 

In this regime, the collisional cross-section is equal to the ge- 
ometrical cross-section augmented by a gravitational focusing 
factor: 

A(i, w) = 7r(i + w)'(l + e) where 9 = ^^Hh^ . (37) 

When the Safronov number 8 is large, this expression simpli- 
fies to 

2 



A(s, Smax) ' 



lirGmraaxS ] 



max max 



(38) 



(J (s) \ Smax , 

In addition, as particles grow larger in size, most of solid ma- 
terial becomes concentrated in fewer and fewer objects, until 
isolation mass is reached (all of the available material is con- 
tained in one object). In this work, I assume that the growing 
protoplanet can indeed accrete all the material available within 
the region of the disk centered on r and of width equal to Ar 
with 

Ar 

= min(-\/A(ii„ax7'Smax),^''H) ■ (39) 

In other words, the total surface density of material available 
for growth (excluding the mass contained in the growing proto- 
planet itself) is 

(40) 



P 27rrAr ' 

Finally, using the expressions derived in §2.4.2 for the parti- 
cle velocity dispersion and for the disk scaleheight, the growth 
of the largest object is found to be governed by the equation 

d-^max ^max'^^max ^TtCs^j., 



df 



'^('^max) 



■/5 



(41) 
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where 



l5- 



(42) 



.Smin/Smu 



and n,nax is reduced to include only the material available for 
growth (see equation ( |40] |) , 



InrAr 



SO that 



dt 



^inax'''-max " 



1.77e- 



nGst 



(43) 



(44) 



lirhp O'(imax) 

2.5.3. Transition size 

The transition between the collisional regime dominated by 
turbulence and the collisional regime dominated by gravita- 
tional interactions is determined by the size for which the es- 
timates of the velocity dispersion are equal, namely when 



<e^(w) >'''^ v'K 



ma> 



(45) 



Note that although this size depends on the surface density and 
temperature of the gas, and therefore on the position within the 
disk, it is typically of the order of a few kilometers. Beyond 
the transition size, the Safronov number is indeed found to be 
much larger than unity, justifying the use of the approximation 
e > 1 in equation (l3BT l. 

2.6. Evolution of the surface density of particles 

The equation of evolution the surface density for each species 
condensed into solid particles is given by Takeuchi, Clarke & 
Lin (2005) for instance, as 

^ + 7a;«+'-s;«p)=o, (46) 

where is the vertically integrated equivalent diffused mass 
flux cause by gas turbulence for each particle species (see equa- 
tion (fT9] l) and Up is the mass-weighted drift velocity of the par- 
ticles resulting from gas drag. 

The radial velocity of a particle of size s was calculated by 
Weidenschilling (1977) and can be written in the notation used 
here as 

u 2TrSt(s) 

"'^'^=4n^StHs)+l-^'''Un^StHs)+l ' ^^'^ 
where 77 is related to the radial pressure gradient in the disk: 

Ih^dlnp 
2 oinr 

Note that the constant 77 reflects the difference between the typ- 
ical orbital gas velocity and the Keplerian velocity at the same 
location in the disk. The mass-weighted average particle veloc- 
ity is then determined by the integral 



2nhr, 



dn 

m(s)up(s) — ds . 
ds 



which integrates to 

Up = Ul ( a/ 27r5fn,ax) " 2r/VK7( V 2vr5fn,ax) , 

where the functions / and J are given by 

/2 

I{x)=^[f,(x) + f2{x)\ and 
\/2 

J{x)= — [-/i(x)+/2(x)] where 



(49) 



(50) 



x^+x 



V2+I 



/i(^)=^ln| 2 , I ' 

2 yx^-xv2+l J 

f2(x) = arctan(xV2+ l) + arctan(x%/2- 1) . 



The functions / and J are shown in Figure [T] Finally, note that 
planetary migration resulting from planet-disk interaction (type 
I or type II migration) is not taken into account here. 



2.7. Sublimation/condensation 

Given the simplistic temperature profile used in this work, a 
simple sublimation/condensation model suffices. The sublima- 
tion and condensation of each chemical species is assumed to 
be instantaneous in time. After each timestep the new surface 
densities in solid and vapor forms are recalculated according to 
the very simple algorithm 



E'(r,f):=S;(r,f) + S;(r,f), 



E'(r,f) 



1 + tanh 



Tj - TJr) 
AT 



E;(r,f) := ^'(r,t)-^(r,t). 



(52) 



(51) 



where 7) is the typical sublimation temperature of the /-th 
species, and AT is taken to be lOK (in practise, the exact value 
of AT only influences the radial extent of the sublimation re- 
gion). 



2.8. Numerical procedure 

The details of the numerical procedure adopted are given in 
Appendix B, for reference. The algorithm constructed follows 
the simple pattern at each timestep, from a given set of initial 
conditions; (i) test whether particles of Size ■^^niax 3^1"^ governed 
by turbulent or gravitational interactions (ii) evolution of the 
particle size though collisions using equations ( l36b or ( l44b ac- 
cordingly (iii) evolution of the gas density (iv) evolution of the 
vapor-phase of each species (v) evolution of the particle phase 
of each species (vi) condensation/sublimation and calculation 
the total surface density of particles. 

The numerical scheme adopted uses a standard split-operator 
techniques, where diffusion terms are integrated using a Crank- 
Nicholson algorithm, the advection terms are integrated using 
an upwind explicit scheme and other nonlinear terms are inte- 
grated using a 2nd order Adams-Bashforth scheme. 

Depending on the spatial accuracy and the number of grain 
species studied, the typical integration time required to evolve 
of a single disk over several Myr varies between a few hours 
and a day on a conventional desktop. 



3. MODEL PARAMETERS AND INITIAL CONDITIONS 

3.1. Model parameters 

The numerical model requires a certain number of input pa- 
rameters, listed in Table 1 ; these are separated between stellar 
parameters, photo-ionizing wind parameters, disk parameters 
and finally grain parameters. Default values for a "fiducial 
model" are also given. 

Table 1 : Fiducial model parameters 
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Fig. 1. — I(x) (solid line) and J(x) (dashed line). As.*: — » 0, 1(x) — » 1 and /(.«:) /3. As x +oo, I(x) an<iJ(x) both tend to V2-ir/4-x. 



Stellar Mass 




IMq 


Stellar Luminosity 






Stellar Radius 




IRq 


Stellar Temperature 


t; 


17b 


Sound speed of ionized gas 


Ci 


10''cm/s 


Amplitude of photo-ionizing flux 




iC^^photons/s 


Turbulent a 




10-^ 


Scaleheight at lAU 




0.0333 


Temperature power law index 




-1/2 


Inner disk radius 


nil 


0.01 AU 


Outer disk radius 


'"out 


2000 AU 


Solid density of grains 


Ps 


1.0 


Sticking efficiency 


e 


10-2 


Separation of protoplanets 


b 


10 



The various values selected for this fiducial model deserve 
comments. The star is chosen to be a solar-type star for ease 
of comparison of the results with the model of Stepinski & 
Valageas (1997) and Ciesla & Cuzzi (2006). Another possi- 
ble choice would have been to select a typical T Tauri star 
(M* = O.SMq, = 4000K, and = 2.5/?©) which was done 
by Dullemond & Dominik (2005). Detailed discussions on 
the values of the parameters associated with the photo-ionizing 
wind can be found in the work of AA07. 

The value of at is selected to be 0.01, which is a reasonable 
upper limit on the value that seems to be favored by numeri- 
cal simulations of MRI turbulence (Fromang & Nelson 2006). 
However, by selecting a constant value of at both in time and 
space, I neglect possible effects of dead-zones (Gammie, 1996) 
which may not exist anyway (see Turner, Sano & Dziourke- 
vitch, 2007) as well as the transition from angular momentum 
transport dominated by gravitational instabilities to angular mo- 
mentum transport dominated by MRI turbulence. The iimer 
disk radius is chosen as a plausible location for the magneto- 
spheric truncation radius (Hartmann, Hewett, & Calvet, 1994) 
while the outer disk radius is chosen at an arbitrarily large dis- 
tance from the central star. 

The solid density of grains ps is an elusive parameter since it 
is quite Ukely to vary strongly with time and with distance from 
the central star, both through repeated compaction events, self- 



gravity (in the case of large objects) and chemical composition. 
Here it is set to unity for simplicity, although this is admittedly 
not very satisfactory. The sticking efficiency is equally difficult 
to constrain a priori, although fascinating computational and 
experimental studies (see the review by Dominik et al. , 2007) 
are beginning to shed light on the subject. Here, I begin by as- 
suming a value of 0.01, and later discuss possible constraints on 
this value from observations of the grain surface density profile 
of disks. 

3.2. Model initial conditions 

The model described in this paper does not take into account 
the evolution of gas induced by self-gravity. It also ignores in- 
fall of mass onto the disk. As a consequence, it is limited to the 
study of disks which are gravitationally stable with negligible 
infall. The "initial" conditions should be thought of as the state 
of the disk after the Class I phase. 

The required initial conditions of the model are: the initial 
surface density of the gas, the initial total surface density of 
heavy elements (both in gas and solid form), the respective pro- 
portion of heavy elements contained in each chemical species, 
and finally the initial maximum size s^ax of the dust particles. 

The initial surface density of the gas is selected to be a trun- 
cated power law (Clarke, Gendrin & Sotomayor 2001) 

S(^,0)=-^e-'-/«°, (53) 



2nRor 

and can therefore be easily characterized by the initial gas disk 
mass Mq = M(0) and the initial disk "radius" Ro- The initial 
total surface density of heavy elements (in both gas and solid 
form) is chosen to be a constant fraction of S(r, 0), with 



Sp(r,0) = ZoS(r,0), 



(54) 



and thus can be characterized by one parameter only, namely 
the initial metallicity fraction Zq. The code is written in a very 
versatile way which allows the user to decide how many sep- 
arate chemical elements to follow. The user needs to input 
the initial mass fraction of each chemical element, as well as 
their sublimation temperature under pressure and density con- 
ditions typical of accretion disks. As a first step, the subUma- 
tion/condensation routine is then run to decide what fraction 
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of the total mass is in solid or in vapor form. The total solid 
particle density is then recalculated accordingly. 

Finally, the initial size of the particles imax('",0) must be cho- 
sen; for simplicity, it is assumed to be constant with imax('", 0) = 
■SmaxO- Although this is clearly an unrealistic initial condition, 
grain growth in the inner disk is so rapid that all "memory" of 
the initial conditions is lost within a few hundred years. On 
the other hand, since growth is negligible in the outer disk, 
Sm!ix(f,t) ~ imaxO there. Hence selecting the value of ^maxO ef- 
fectively determines the timescale for the evolution of solids in 
the disk (see ^6. 2b . While the fiducial model considers imaxO 
to be equal to the maximum plausible particle size in the MRN 
size-distribution function for the ISM, one could also imagine 
grains to grow even in the core-collapse phase. Suttner & Yorke 
(2001) found that grains could achieve sizes up to 10/j,m post- 
collapse, and so I will consider cases with varying initial con- 
ditions for imaxO in addition to the fiducial model (see ^5.31 ). 

Table 2 summarizes the initial condition input parameters, 
and gives typical values for a fiducial run. 

Table 2: Fiducial model initial conditions. 



Initial disk mass 


Mo 


0.05M,, 


Initial disk radius 


Ro 


30 AU 


Initial metallicity 


Zo 


10"^ 


Number of species 


maxtype 


3 


Initial i^^x 


■^maxO 


1/im 



The initial chemical composition of the dust, in the fiducial 
model, is taken to be the following: 45% "Ices" and other 
volatile materials (with sublimation temperature Tic = 170K), 
35% refractory material (with sublimation temperature T^i = 
470K) and 20% finally iron-based material (with sublimation 
temperature Tpe = 1300K). The solid composition and subli- 
mation temperatures are adapted from Table 2 and Table 3 
of Pollack et al. (1994) to account for a reduced number of 
species. 

The fiducial initial model (after condensation/sublimation of 
the relevant species) is presented in Figure |2] 

3.3. Model tests 

The numerical algorithm was tested against the results of 
AA07 for the evolution of the gas and grains by using their 
initial conditions, switching off grain growth, sublimation and 
condensation, and by replacing equation ( |50l l for the drift veloc- 
ity with equation ( |47] i. Both gas and grain evolution are found 
to be in perfect agreement, as required. 

4. OVERVIEW OF RESULTS IN THE FIDUCIAL MODEL 

The fiducial model presented in ^ was integrated forward 
in time until complete dispersal of the gas. Figure [3] shows the 
evolution of the surface density of the gas, the total solid sur- 
face density as well as that of the three species considered. Fig- 
ures]?^ andU^ show the evolution of the particle size and total 
metallicity as a function of radius and time. Finally, Figure |5] 
shows the evolution in time of the total mass of gas and dust in 
the disk. 

4.1. Evolution of the gas surface density 

The characteristic evolution of S(r, f) under this particular 
photo-ionizing wind model has been extensively studied by 



Alexander, Clarke, & Pringle, (2006a and 2006b) (see also 
Clarke, Gendrin & Sotomayor, 2001). It can be seen in Fig- 
ure [3] as a dotted line, and in more detail in Figure |6] While the 
mass flux from photo-evaporation is negligible compared with 
the mass flux from viscous accretion/spreading, the disk under- 
goes a long period of near self-similar evolution. When both 
fluxes become comparable a depression appears in I](r, f ) and a 
gap eventually forms, here at radius rgap = 0.9 AU, at f = 7Myr 
Within a few thousand years, most of the gas in the inner disk 
has been accreted onto the central star, while the radius of the 
hole begins to expand as a result of direct photo-evaporation. 
At f = 7.12Myr, the hole radius has retreated to 200 AU, and 
finally beyond 500 AU after t = 7.19Myr 

While the evolution of the gas is (in this model) independent 
of the evolution of solids, particle growth and particle migration 
are nonlinearly strongly coupled. 

4.2. Particle growth 

The evolution of the maximum particle size imax('", is 
shown in Figure|4^ both for very early times and at later times. 
Grain growth is extremely rapid in the inner disk regions in the 
early stages of disk evolution, in particular near sublimation 
lines. Within just 100,000 yr, a characteristic shape to the curve 
SmsLx(t';t) appears, which contains three different regions: (I) in 
the innermost disk region (r smaller than a fraction of 1 AU), 
a slightly tilted plateau corresponding to the particles having 
reached isolation mass; (II) a power-law region (from a frac- 
tion of 1 AU to about lOOAU); (III) a region where grains have 
undergone negligible growth. Superimposed on this character- 
istic shape are a set of peaks corresponding to the successive 
sublimation lines. The transition between region I and region 
II is easily identified as the transition between the gravitational 
regime and the turbulent regime; its steepness confirms that as 
soon as gravitational focusing becomes effective, the collision 
rate increases and particles rapidly reach isolation mass. The 
transition between region II and region III can also be easily 
identified as the region where the growth timescale of particles 
of size imaxO bccomcs comparable with the age of the disk. 

Once established (after the first Myr), the global shape of the 
curve imax('",0 varies little with time (see solid lines), although 
particles within the sublimation region continue growing, and 
the three regions slowly expand outward. 

4.3. Solid density and chemical composition 

The evolution of the solid density is shown in Figure[3] Small 
particles well-coupled with the gas (5fniax ^1) closely follow 
its inward or outward motion depending on the radial position 
considered. As a result, during the initial viscous spreading of 
the disk (within the first Myr) a significant proportion of the 
mass in solids is transported outward with the gas creating a 
large reservoir of small dust grains at large radii (r > lOOAU). 
Meanwhile, particles in the inner regions of the disk (r < 10 
AU) rapidly grow and begin to drift towards the central star dif- 
ferentially from the gas, which results in local changes in the 
metallicity Sp/S. Figure |4]3 shows the evolution of metallicity 
in more detail, and reveals that the inner and intermediate disk 
go through an initial phase of strong depletion in heavy ele- 
ments. Later, the global evolution of the surface density of par- 
ticles is controlled by the mass flux incoming from large radii. 
The observed increase in the metallicity is essentially related to 
the decrease in E through photo-evaporation. 

In addition to this global trend strong surface density peaks 
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Fig. 2. — Initial dust and gas surface densities in the fiducial disk model. The dotted line corresponds to the molecular gas and the solid line to the total surface 
density of solids. The three species considered are: the volatile material (dot-dot-dot-dash line), refractory material (dot-dash line) and the iron-rich material (dashed 
line). 



can be observed near the successive sublimation lines. These 
are presumably caused by the differential drift of the solid and 
vapor form of each chemical species (Stepinski & Valageas 
1997, Ciesla & Cuzzi, 2006). The peaks consistently stand 
roughly one to two orders of magnitude above the smoother 
"background" surface density profile, but do not appear to grow 
independently of it beyond the first 100,000 yrs. As the gas den- 
sity decreases, the local metallicity near the sublimation lines 
steadily grows. By 4Myr, the sublimation line for refractory 
materials has equal content in gas and solids, suggesting the 
possibility of local onset of gravitational instability of solids 
(which is not modeled here). After 6Myr, the two remaining 
sublimation lines (for the volatile and iron-rich material) also 
pass the same threshold. Interestingly, Figure |3]reveals that the 
very rapid growth of material near each sublimation line traps 
a variety of grain species into the growing bodies, so that the 
strong enhancement in the surface density of icy bodies near the 
volatile sublimation line is accompanied by an enhancement in 
the surface density of refractory materials and iron-rich mate- 
rials. The same phenomenon is observed near the sublimation 
line for refractory materials. 

4.4. Disk masses 

The evolution of the total mass in gas and solids in the disk 
is shown in Figure|5]as dotted and solid lines respectively. Also 
shown are the total amount of solids found within 20 AU and 
outside of 20AU. At f = 0, the solid mass is equally distributed 
between the inner (< 20 AU) and the outer disk (> 20 AU); 
Within a short time (of the order of 100,000yr), most of the 
mass in the inner disk accretes onto the central star, while the 
mass contained in the outer disk remains at a constant fraction 
of the disk mass in gas. Beyond this point, the mass in the inner 
disk is controlled by the flux of material drifting in from the 
outer disk. 

When the gap opens (at 7Myr), the total mass of gas drops 
precipitously (within about 200,000 yr) while the total mass of 
solids remains constant. The total content of solids left in the 
disk after complete photo-evaporation of the gas is about 1.3 
X 1O~^M0, or in other words about 4 Earth masses. Only 20% 
of this amount is located within the inner 20AU of the disk. 



while the remaining 80% are swept out to the outer disk. 

5. MATHEMATICAL INTERPRETATION OF THE RESULTS 

In order to gain more insight into the numerical results for 
the fiducial model, it is useful to characterize and when pos- 
sible quantify some of the generic behaviors observed in the 
solutions. 

In this section 1 present both existing and new analytical re- 
sults on the evolution of gas and solids in viscously evolving 
disks. While the complexity of the system clearly precludes the 
existence of a closed-form fully analytical solution, there are 
certain limits where analytical efforts pay off. By comparing 
the analytical estimates derived with the exact outcomes of the 
numerical algorithm, 1 am able to test the numerical results on 
a systematic basis, and at the same time obtain strict constraints 
on the regimes of validity of the analytical solutions. 

Given that the evolution of the gas is more-or-less indepen- 
dent of the evolution of solids, much progress has already been 
done in describing it analytically. These are presented in ^5.11 
New results on the evolution of solids are presented in ^5.2| and 

£3] 

5.1. Evolution of the gas 

The evolution of the gas density is shown in more detail in 
Figure |6] Lynden-Bell & Pringle (1974) (see also Hartmann 
et al. 1998) showed that provided (i) the mass accretion rate 
due to viscous transport is much larger than the wind photo- 
evaporation rate and (ii) the disk is allowed to spread to infinity, 
then there exist a simple self-similar solution for S(r,f): 

E(,,0=:^r-V2expf 



T=Ul 

TV 



with 



(55) 



where the viscous spreading time of the disk is Tv = R^/1>Vi{Rq). 
In this solution, the gas velocity is equal to 

'--—) 

showing that m < for r < r^{t) (viscous accretion) and m > 
for r > ry{t) (viscous spreading), the critical radius being r^{t) = 



(56) 
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RqT /2. Figure|6]compares the self-similar solution compares 
with the true numerical solution, and reveals excellent agree- 
ment at early times, gradually deteriorating in the inner disk 
as time evolves and photo-evaporation begins to dominate the 
gas dynamics. In the outer disk the agreement remains glob- 
ally much better (since the photo-evaporation rate is very low 
at large radii) and deteriorates only slightly (by a factor of no 
more than a few) as expected when the critical radius reaches 
the outer boundary. Note that a much better approximation to 
I](r, t) including the effects of photo-evaporation has been ob- 
tained by Ruden (2004). 

In the early self-similar phase, and within r^, the mass ac- 
cretion rate M{r,t) = lirruYi is roughly constant with radius. 
The total gas disk mass decays as M{t) = MqT'^I^ (neglect- 
ing the effects of the outer boundary condition), so that the gas 
accretion timescale increases linearly with the reduced time: 
M/\M\~2Tt^. 



When the photo-evaporation rate becomes comparable with 
the accretion rate, a gap opens in the disk. HoUenbach et al. 
(1994) argued that the gap opening radius rgap is located close 
to the gravitational radius rg = GM*/c„ while Liffman (2003) 
and Font et al. (2004) revised this estimate to be a fraction of 
r„. Since r„ scales linearly with stellar mass so does rgap; in 
what follows, I adopt 

The time at which the gap opens can be estimated by equating 
the wind mass loss rate Mw to the viscous accretion rate in the 
self-similar solution (Clarke, Gendrin & Sotomayor, 2001); this 
yields (provided Tgap ^ Tv) 



2tvM„ 



2/3 



(58) 



The wind mass loss-rate prior to gap opening for the fiducial 
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Fig. 4. — Left: Evolution of the maximum particle size at selected times. From bottom to top, ( = 0, 10** (dotted line), lO' (dash Hne), 10* (dot-dash line), 2 X 10* 
(dot-dot-dot-dash Hne), 4 X 10* (long-dash line) and 6 X 10* yr (solid Hne). Note the strong growth peaks located near the respective subHmation lines, the plateau 
for r < 0.1 AU where particles have reached isolation mass and the region of negligible growth for r > lOOAU. Right: Metallicity fraction at the same selected times 
as in the left-hand- side figure. Note the strong initial reduction caused by the rapid inward drift of the particles, followed by gradual growth. The latter is caused by 
the reduction in E rather than by an increase in Sp. 
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Fig. 5. — Total disk mass in the fiducial model. The dotted line shows the integrated gas mass and the solid Hne shows the integrated solid mass. To illustrate the 
rapid loss of solids in the inner disk, the total disk mass contained in r < 20AU is shown in the dashed line, while the rest is shown in the dot-dash line. 



model was calculated by AA07 to be 



Mw~4x 10"'"MQ/yr. 



(59) 



Table 3 compares the estimate from equation (l58T l with the 
outcome of numerical simulations with varying Rq and Mq, 
showing that it is indeed a very good estimate for the gap for- 
mation timescale except for the fiducial model. This is because 
the actual viscous mass accretion rate in the numerical solu- 
tion of the fiducial model deviates from the simple estimate of 
M()T~^^^ /It^ at large times, when the disk spreads all the way 
out to the outer edge of the numerical meshQf 

Table 3; Gap opening time, actual and predicted 



Ro (AU) 


Mo/Mq 


r^T (Myr) 


(Myr) 


30 


0.05 


7.01 


7.79 


10 


0.05 


5.35 


5.40 


5 


0.05 


4.36 


4.29 


10 


0.01 


1.96 


1.84 



After the gap formation, supply of material from the outer 
disk is shut off and the inner disk rapidly clears of all gas. The 
gas clearing time can be estimated from the viscous timescale 
at the gap formation radius: 



'^clear ~ '^v 



gap 



(60) 



In the fiducial model, the gas clearing time is of the order of 
4,000 yr only, and can be considered to be near-instantaneous. 
This is indeed seen in the simulations (see Figure[3]). 

* To get a better estimate for the gap formation time in this case, one should simply not neglect the outer boundary term in the calculation of the total mass of the disk, 
see Hartmarm et al. 1998. 
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Fig. 6. — Gas surface density in the fiducial model (from the top to the bottom curve) at f = 0, 1 , 2, 3, 4, 5, 6 and 7 Myr The solid lines show the exact numerical 
solution, while the dotted lines show the analytical estimate provided by the self-similar solution i55\ at the same times. 



Direct photo-ionization of the gas at the hole edge results in 
a sharp change in the gas mass loss rate (see Appendix A). The 
evolution of the size of the hole can be derived from the work 
of Alexander, Clarke & Pringle (2006b) to be 

2 



?"hole(f) = rg. 



2t, 



+ 1 



clear 



(61) 



which is again a good estimate of the hole radius derived from 
the numerical solution (see ^4.11 ) for t > Tg^p. 

5.2. Grain growth 

For simplicity, in this section I focus on deriving analytical 
estimates for grain growth in disks where sublimation and con- 
densation are neglected. For ease of comparison with the nu- 
merical model, an additional run was performed using the fidu- 
cial model but with no sublimation or condensation of mate- 
rial, yielding the solution shown on Figure |2l After the initial 
rapid growth period (for t > 10^ yrs), one can note very clearly 
the three regions of interest described earlier: region III where 
Smax(r,t) ~ imaxo; region 11 where .Wx('",f) appears to follow 
roughly a power law in r, and region I where the particles have 
reached isolation mass and stopped growing. In this figure the 
peaks associated with sublimation lines are naturally absent. 

Despite the apparently simple behavior of the numerical so- 
lutions, modeling this evolution is a complex problem: the par- 
ticle growth rate depends on Y,p(r,t), which is regulated by the 
drift velocity of the particles, which depends on the particle size 
Smix{f,t)- The implied nonlinearities preclude the existence of 
analytical solutions in most cases. 

However, there exist a limit in which insight can be gained 
from simple models: when the particle drift time is much longer 
than the particle growth time, one can expect the metallicity to 
remain close to its initial value, namely Sp/S = Zq. This hap- 
pens in the very early stages of disk evolution (the limit of small 
f ), as well as in the outer regions of the disk (region III) at all 
times. 

5.2.1. Growth timescales 

To interpret the numerical solutions I consider the growth 
regime dominated by turbulent encounters, where particles of 



size imax follow the growth law given by equation ( [36] l. The 
growth timescale of the particles is given by 

-I 

(62) 



1 dSn 



In the limit where St„ 



V at 

while in the limit where 5fmax ^ 1 



dt 

^ 1, the growth timescale is equal to 



Td 



(63) 



(64) 



Note that the second expression is independent of a^. This is re- 
lated to the fact that for larger particles, coupling with the turbu- 
lent eddies is weakened, which reduces their relative velocities 
but also increases particle concentration through sedimentation. 

Figure |8] shows (as lines) the exact analytical expression for 
the growth timescale obtained by combining equation ( l62b with 
equation ( |36] |, and considering Sp/S fixed and equal to the ini- 
tial metallicity Zq = 0.01 . Growth timescales of particles at radii 
0.1, 1, 10 and 100 AU are shown. The power laws seen in ei- 
ther limits are well-approximated by equations ( |64l ) and ( |63] l, 
and the flattening of the four curves corresponds to the transi- 
tion 5finax ^ 1 to 5fmax ^ 1 ■ This figure is particularly useful 
for reading directly the maximum size of particles achievable 
at a given age in the disk should the surface densities of dust 
and gas indeed remain constant over that period. For exam- 
ple by looking at the intersection of a horizontal line at 10"^yr 
with each curve, one deduces that in a lO'^yr-old disk there will 
be no significant growth beyond 20 AU, 10-20/im-size parti- 
cles at lOAU, m-size objects at 1 AU and finally, gravitation- 
ally dominated growth (towards isolation mass) below 0.1 AU. 
Since Tg oc I/Zqe, the growth timescale and maximum particle 
sizes for disks with other values of the metallicity Zq or stick- 
ing efficiency e can be read by translating the curves up or down 
accordingly. 

In reality, the metallicity Z = Sp/S is of course not con- 
stant. Figure [8] also shows (as symbols) the actual particle size 
achieved in the fiducial disk with no sublimation/condensation 
at the same selected radii. The lines and the symbols follow 
each other reasonably well in the expected limits (large radii or 
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Fig. 7. — Maximum grain size iji^axC'^:^) at 7 — 0, 1. 10, 100, 1000, 10^, 10^ and finally 10*^ yrs (solid lines, from bottom to top) obtained from the numerical 
simulation of the fiducial model with no sublimation/condensation. The dotted lines show the analytical estimates of equation (65) at the same times. The dashed 
line marks the transition between the turbulent and gravitational regime axt = IMyr 
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Fig. 8. — Growth timescale for particles of size .v for the initial conditions of the fiducial disk described in j|3.2l with no sublimation/condensation; the curves 
show the growth timescale at r = 0. lAU (solid line), r = lAU (dotted line), r = lOAU (dashed line) and finally r = lOOAU (dot-dash line). The symbols show the 

true A'max 

(r,0 at the same radii (0.1 AU: diamonds; 1 AU: stars; 10 AU: triangles; 100 AU: squares). 



short times) up to the point where 5fmax — 1 at which point the 
simple analytical estimate systematically breaks dowrQ. 



5.2.2. Particle size 

Given the good fit found for particles with 5fn,ax <C 1 , 1 now 
approximate s^^,^(r,t) by the value of the grain size for which 
the growth timescale equals the age of the disk: setting Tg = f in 
equation yields ( |63] ) 




Ps \ lyear y V / V-"'^© 

'(65) 

If one assumes as before that Sp/S ~ Zq then one can get a 
rough estimate of s^^^ir, t) by combining equations ( |65] ) and 
( |55] ), and shown in Figure|7]as dotted lines. 



As expected, the estimate for SmaxC?",?) is in fairly good a 
agreement with the numerical results for small times; it cor- 
rectly predicts the power law structure of the whole intermedi- 
ate region (for s < 1km, roughly) at early times (t < lO'^yr), but 
not so at later times (where the analytically predicted power law 
is too steep compared with the numerical results). This is again 
related to the fact that the surface density of particles becomes 
significantly depleted at later times. 

The estimate for imax('", also correctly predicts the transi- 
tion between regions 11 and 111 of the disk. Since particle growth 
is fundamental to our understanding of disk SEDs, 1 now give 
an analytical estimate for this transition radius: for early times 
(for t < Tv), rjjj(t) is given by 

1 1/4 / . X 1/2 / X 1/4 



'm(0 = 



Z^e^Ui Rq 



St 



maxO 



lAU 



lyear 







AU, (66) 



^ The reason why Sfmax ~ 1 is equivalent to the point where Z begins to differ significantly from Zq is related to the fact most of the mass is contained in the particles 
of size i'max^ which also happen to be the particles with the highest inward drift velocity. 
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where 



St 



maxO " 



(67) 



/IttjMq 

is the Stokes number at f = and r = Roof particles of size Sm^xO 
For later times (t > t^) 



'■m(t) = 



5fmax0 lAU 



1/4 



\ lyear / 



3/8 



/ t 



\ lyear 



1/8 



Ml. 

Mq 



1/4 



AU 



(68) 



5.2.3. Gravitational regime 



The transition from the turbulent regime to the gravitational 
regime (region II to region I) is easily understood by consider- 
ing equation (05]). The transition size S^-dx is given by 



W(r,f) = 7.1( — 



10/51 /^^^ 8/17 



I](r,f) 



7/17 



lOOOg/cm^ 



km. 



(69) 

The curve for Ssmax('", at f =lMyr is shown on Figure |2l and 
correctly marks the transition between the turbulent and gravi- 
tational regime at the time considered. As time progresses and 
I](r, f ) decreases so does the transition size. 

5.3. Evolution of the solid mass fraction prior to gap opening 

The evolution of the solid mass fraction is governed by par- 
ticle diffusion and drift. The analytical prescription used to de- 
scribe the particle size distribution function is particularly use- 
ful since it can easily be integrated to yield the bulk motion 
properties, as seen in equations (|20] | and (l50l l. Given the asymp- 
totic behavior of the functions I{x) and J{x), Up is roughly equal 
to 

Mp = u-lriva for Stm^x < 1 , 



gradually decouple from the gas implying that the reservoir be- 
gins to "leak". Eventually, even the smallest particles decouple 
from the gas, and all the solids come rushing inward. The evo- 
lution of the particle velocity can be seen in Figure|9] Note the 
existence of the reservoir (a region of significant radial extent 
with Up > 0) for all particle sizes at early time. As time pro- 
gresses, the reservoir begins to "leak" as the particles gradually 
decouple from the gas, until a a critical point where Up becomes 
negative for all radii. The phenomenon depends strongly on the 
initial size of the particles Sm-dxO- three simulations are shown 
in which the maximum particle size is respectively l/im, 3/im 
and 10/im. The timescale tp for the release of the particle reser- 
voir is clearly much shorter for the larger particles (0.84Myr for 
10/xm-size particles instead of 2.33Myr for /im-size particles). 

In fact, this timescale can easily be estimated by solving si- 
multaneously the equations Mp = and dup/dr = 0. Assuming 
that I](r, f ) is equal to the self-similar solution, and that the par- 
ticles in that regime satisfy 5fmax ^ 1 (which was checked nu- 
merically). Up ~ M-47r775fmax/3. It follows that 



+ 1 = 7p = 



TdiRo) 1 



-,2/5 



leeir^rjiRo) 



St, 



maxO 



(71) 



where rd(/?o) is the dynamical time at Rq, 77(/?o) is obtained by 
applying equation (l48T l at r = Rq, S'fmaxO is the Stokes number 
of particles of size imaxo at Rq (see equation (l67Ti). Note that 
in order to derive this expression, I have also made explicit use 
of the fact that q = -l /2. The quality of this estimate is found 
to be excellent given the approximations made (see Figure [TOli. 
and small discrepancies are attributed to the fact that S(r,f) is 
not exactly equal to the self-similar solution, and that Mp has 
been approximated by its Taylor expansion for small St max- In 
the same analytical calculation, the radius of the reservoir at re- 
lease is found to be = R^Tj,; checking this solution against the 
numerical runs also reveals excellent agreement. 



Mp = (m-2?7Vk) 



for Stmdx > 1 



(70) 



so that, as expected, the bulk radial velocity of the particles 
is close to that of the gas for small particles, and tends to 
when imax grows. Note that Up oc 1 /\^St^ instead of 1/St(s), 
which accounts for the fact that even though particles of size 
imax may be largely decoupled from the gas, a non-negligible 
mass fraction is contained in rapidly drifting intermediate-size 
particles. This is the main difference between this model and 
a single-size particle model; it accounts for the fact that even 
when max ^ 1 , a significant fraction of the coUisional encoun- 
ters are destructive and result in the erosion of the larger bodies 
into smaller rapidly drifting particles with St{s) ~ 1. 

5.3.1. Reservoir of small grains at large radii 

The sign of Mp determines whether grains are transported in- 
ward or outward. As expected from equation ( |56] | the gas ve- 
locity changes sign at ry{t). This critical radius grows linearly 
with T (defined in equation (|55]l), and therefore sweeps out- 
ward roughly on the viscous timescale t^. As a result for r < r^ 
the particle velocity is necessarily negative, while for r > r^ 
particles can be entrained outward provided they are strongly 
coupled with the gas. Since r^ is typically much larger that r{{j, 
all the particles outward of r^ are particles of size s^axu- Com- 
bining these facts implies that there exists a reservoir of small 
particles at large radii, slowly eroded by the outward motion of 
r^it). In addition, as the gas density drops, the small particles 



5.3.2. Evolution of the total mass of the particle disk 

Since particles in the inner disk can achieve significant sizes, 
the drift timescale of grains within rv(f ) is usually much smaller 
than the viscous accretion timescale and/or the age of the disk. 
This can be readily seen in Figure|9l 

This very simple fact has important consequences: it implies 
that the distribution of solids in the inner disk is uniquely con- 
trolled by the mass flux leaking from the reservoir. One way 
to see this is to look at the particle disk evolution timescale 
Mp/|Mp| obtained by numerical integration of the fiducial 
model. Figure [TT] shows results for the three different initial 
particle sizes considered in the previous section (as solid lines). 
The linear increase for early times (T ^Tp) mirrors the gas evo- 
lution timescale, as expected from the tight coupling between 
the reservoir particles and gas: 



Mpit) = Z^it) = -^T-'l\ 



(72) 



so that Mp/\Mp\ ~ M/M = 2t^T . As T exceeds Tp, the hnear 
increase saturates then rapidly turns over, as expected from the 
release of the reservoir particles. 

Constructing an exact analytical model governing the evolu- 
tion of the particle disk after Tp from first principles turned out 
to be rather difficult. However, it is possible to gain insight into 
the problem by inspecting the results of the numerical simula- 
tions first. At later time, one can expect the dynamics of the 
particle disk to depend on the reservoir release timescale Tp. I 
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Fig. 9. — Radial velocity u of the gas (soUd line) and mass-weighted average radial velocity Mp of particles of maximum size Smaxo = 1m™ (dotted line), 3/im 
(dot-dashed line) and lO^xm (dashed Une). Note how the slow erosion and eventual release of the reservoir (i.e. of the spatial region with up > 0) depends on imaxo- 
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Fig. 10. — Non-dimensional timescale Tp for the release of the particle reservoir as a function of initial particle size ^maxo- The analytical estimate (solid line) is 
compared with hand-checked values of Tp for three numerical simulations of the fiducial model with different initial particle sizes (stars). On the same plot is shown 
the numerically determined radius of the reservoir at release rp in AU (diamonds) as well as the corresponding analytical estimate RoTp (dashed line). 



seek a functional form of the kind 

Mr, ' \r ' 



(73) 
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with 



1 + a2x" 



(74) 



within the disk is controlled by the mass flux coming from the 
reservoir one can expect that 



(with ai, 02, fl3 > 0) which satisfies the requirement that f(x) 
1 as X ^ 0, and f(x) nearly exponentially as x oo). A 
fairly good (but clearly not perfect) fit for all three curves is 
empirically found to have ai = 0.2, a2 = 0.25 and = 3.3, im- 
plying 



dMp Mp / T 
~ ^-exp 0.2 — 




(75) 



The fitting curves, for each of the three initial particle sizes cho- 
sen, are also shown in Figure [TTh. The initial linear rise as well 
as the maximum are very well represented, while the fit at later 
times (for T > 7p), in particular for the smaller particle sizes, 
is slightly poorefl Integrating equation ( fTSl ) yields an estimate 
for the total particle disk mass as a function of time, which 
is compared with the results of the numerical simulations in 
Figure [TTb . The approximate solutions follow the trend of the 
exact solutions well, with some small acceptable systematic de- 
viations at early times (see below). In particular, it reproduces 
well the very rapid decrease in the particle disk mass as the 
reservoir is finally released (T > Tp). 

A very important consequence of the "leaky reservoir" model 
is that the total disk mass is reasonably independent of the phys- 
ical phenomena taking place in the inner disk (provided the bulk 
drift timescale of the particles is smaller than the viscous accre- 
tion timescale). This implies that the evolution of the total disk 
mass depends very weakly on the particle growth rate (and in 
particular of the sticking efficiency e), and of sublimation or 
condensation fronts. This can actually be seen in Figure [TT] 
which shows the numerical solution for the fiducial model with 
sublimation/condensation in addition to the three curves dis- 
cussed earlier. The disk mass in the fiducial model is practi- 
cally indistinguishable from that of the disk without sublima- 
tion/condensation (for the micron-size particles). 

Finally, note that the proposed evolution equation for Mp 
breaks down if the bulk drift timescale of the particles is com- 
parable to or larger than the viscous accretion timescale or the 
age of the disk (whichever is smaller). For instance, if the par- 
ticles remain small at all times, then they naturally follow the 
evolution of the gas at all times (which explains the results of 
AA07). As an other example, one can see in Figure [TT] when 
T I that there is a very small difference between the fidu- 
cial model with sublimation/condensation lines and the model 
without. This arises because at early times, the drift timescale 
of the particles is much smaller than the age of the disk and 
the solid mass content has not yet had time to equilibrate. As 
a result, most of the solid mass in the inner disk rapidly drifts 
toward the central star. In the absence of sublimation lines all 
of this mass is lost, whereas a significant fraction of it can get 
trapped by the sublimation lines if subhmation/condensation is 
taken into account (see ^5.51 ). 



5.3.3. Evolution of the particle surface density 

Having characterized the global evolution of the total disk 
mass in solids (which was shown to depend only on the initial 
grain size ^maxO and on the viscous diffusion time Tv), one could 
hope to describe the evolution of the surface density of grains as 
well. As mentioned earlier, when the surface density of grains 

^ This is again attributed to the numerical boundary effects plaguing the gas disk 



27rrMp(r,f)Sp(r,f)-Mp(f), 



(76) 



where Mp(f) is given by equation dTSl ). This approximation 
turns out to be quite good (except outside of r^it) of course). 
Unfortunately, the problem lies elsewhere: even though Mp(f) is 
known, it is particularly difficult to estimate Mp(r,f) since it de- 
pends on J'max('", t) and S](r. t): the analytical estimate of J'max('", 
given in equation ( |65] ) is unfortunately not valid in regions I and 
most of region II beyond a few times 10"* yr and the self-similar 
solution for I](r,f) is also invalid for f > Ty in the same regions. 
Using these estimates despite their poor quality yield predic- 
tions that are off by up to an order of magnitude (see Appendix 
C). Note that if all that is needed is a "quick and dirty" order 
of magnitude estimate then the procedure described in the Ap- 
pendix could be considered satisfactory. In particular, it does 
reproduce well the particle surface density dip observed in the 
intermediate disk regions, a feature which could be used to- 
gether with spatially resolved disk observations to constrain the 
value of the sticking efficiency e. 

5.4. Evolution of the solids posterior to gap opening 

The study of the evolution of solids after the gap opening 
phase was the main purpose of the work of AA07. When 
growth is ignored, AA07 showed that the evolution of the sur- 
face density of solids follows that of the gas (for small parti- 
cles). 

After the formation of the gap, both gas and solids within 
quickly accrete into the central star on the clearing timescale, 
leaving a hole clear of both dust and gas. Near the edge of 
the hole thus formed, an inversion in radial pressure gradient 
causes the particles to drift outward instead of inward. As the 
hole grows all of the grains outside of r^okit) are slowly shep- 
herded outward with it. As a result of these processes, most of 
the solid mass at the time of gap opening is retained in the disk 
but moved to larger and larger radii. 

When grain growth is taken into account on the other hand, 
a different picture emerges. The particles within the initial gap 
radius have typically grown to embryo sizes, so that their Stokes 
number is well above unity by the time the gap opens. The re- 
duction in the surface density of the gas does not affect their 
drift velocity much (which is exactly the opposite case to the 
AA07 model) and they stay in place while the gas within the 
gap accretes. The "hole" thus created still contains a significant 
amount of solid material. 

As direct photo-evaporation takes over, the hole in the gas 
widens as expected but the corresponding reverse radial pres- 
sure gradient has very little effect on the large particles. As a 
result, one observes a significant amount of solid material re- 
maining in the inner disk all the way out to about 10 AU (in 
the case of the fiducial model). Eventually, the edge of the hole 
retreats out to regions where the particle Stokes number is of 
order unity, at which point the clearing begins as in the AA07 
model. As akeady suggested by AA07, all particles smaller 
than 1-lOcm are entrained to large radii, while all particles 
larger than this particular size remain behind. 

The prediction for the particle size given in equation (l65b can 
be used to estimate the radius outside of which particles are en- 
trained, by setting the age of the disk to be f = Tg^p: let rinner the 

evolution for ; > 6Myr, which introduce some non-self-similar effects in the solution. 
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Fig. 1 1 . — Left figure: Particle disk evolution timescale Mp/\Mp\. The solid lines are the outcomes of the numerical simulations for the fiducial model with initial 
particle sizes (from top to bottom) l^m, 3fim and lO^m respectively. The empirical analytical fits to these formula are also shown, with the dotted line for the Ifim 
case, the dot-dash line for the 3/im case and the dashed line for the lO^m case. Right Figure: Total mass in solids in the disk as a function of the nondimensional 
time T. The three solid lines are the outcomes of the numerical simulations for the fiducial disk model with initial particle sizes (from right to left) Ifim, 3fim and 
lO^im respectively. Approximates to these numerical results are obtained by seeking the solution to equation il5i . and also shown here (linestyles are the same as 
in the left-side figure). 



innermost extent of the cleared region, so that 



1/3 



AU 



(77) 



lyvj » \Mq 

which yields a prediction of Hnner =16 AU for the fiducial 
model. This value is larger than that observed in Figure [3] by 
a factor of about 2 ~ 10'/"', which can be attributed to the fact 
that the estimated particle size at this point is a factor of 10 too 
large compared with the true numerical simulations (see Figure 
|4]l. However, the estimate in equation ( iTTb can be thought of as 
a solid upper limit for the radial extent of the remaining solid 
material in the inner regions of the disk after the clearing of the 
gas. 

Material outside of Hnner is shepherded outward with the re- 
treating hole. As more and more material is swept and en- 
trained outward, a strongly localized surface density peak ap- 
pears, which eventually grows to be a large as the local surface 
density of the gas. When this happens, two situations could 
occur within similar outcomes: either gravitational instabilities 
set in, resulting in the in-situ formation of a planet, or the fric- 
tional drag between the particles and the gas begins to influence 
the evolution of the gas itself, and the particles stop moving out- 
ward (neither effects are included in the numerical model, and 
therefore cannot be seen in the simulations). This effect is dif- 
ficult to quantify in the general case, since neither Sp nor E are 
known analytically at this stage of the disk evolution. In the 
fiducial model show in Figure [3] this occurs at about 200AU. 
This process could lead to the systematic formation of local- 
ized debris rings reasonably far from the central star (from a 
few tens to a few hundreds of AU), without any need for prior 
planet formation or other clumping mechanism. Such debris 
rings are commonly observed, or inferred from the dust dus- 
tribution (see for instance Schneider et al. 2006 for direct de- 
tections, and Strubbe & Chiang 2006 for inference from spatial 
dust distribution). 

As in the AA07 model, most of the solid content present in 
the disk at f = Tgap remains in the disk after complete clearing 
of the gas. The fraction of solid material moved to large radii 
depends on the details of the surface density distribution of par- 



ticles at Tgap (which is not well-known a priori), but is typically 
of the order of 80%-90% of the total mass of solids; the remain- 
ing fraction can be found in the inner disk. 

5.5. Effects of sublimation and condensation 

The role of sublimation lines on the local accumulation and 
growth of particles has already been shown and discussed by 
others (Stepinski & Valageas, 1997; Ciesla & Cuzzi, 2006). 
Roughly speaking, the idea is that particles composed mainly 
of a given chemical species drift inward until they reach the 
sublimation line where they are transformed into vapor The 
vapor diffuses inward much more slowly than the rate of mi- 
gration of the incoming solids, leading to a large enhancement 
of the local metallicity. Through turbulent diffusion, a fraction 
of the vapor content actually finds it way back through the sub- 
limation line and recondenses. The typical width of the region 
where this effect dominates can be evaluated from a local diffu- 
sive lengthscale, and naturally scales as h{r). The exact amount 
of material accumulating in the region is more difficult to es- 
timate a priori, since it depends on the difference between the 
solid mass flux into the sublimation line and the vapor mass flux 
out of the sublimation line. However, for most of the lifetime of 
the disk the flux of both solid and gaseous material is controlled 
by the reservoir at large radii so that the solid mass flux is equal 
to ZqM, which is also roughly equal to the mass flux of va- 
por out of the sublimation zone. This explains why the relative 
strength of the surface density peaks compared with the back- 
ground curve fails to grow with time after the initial adjustment 
period, which would necessarily occur otherwise. 

The conclusion is that any local surface density enhance- 
ment, and associated localized peak in the particle size is cre- 
ated very early in the lifetime of the disk (this is indeed seen in 
Figure ^) - but it is only later, when the gas surface density 
decays and the metallicity increases, that gravitational instabil- 
ities in the particle layer could set in to trigger the planetary 
formation process. This picture, should it be correct, also im- 
plies that the relative heights of the peaks determines a strict 
sequence of "alarm clocks" on planetary formation timescales. 

Note that only three species have been selected in the fiducial 
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model. Clearly, there will be as many surface density and par- 
ticle size peaks as the number of separate sublimation temper- 
atures. Also, in this particular model the background tempera- 
ture of the disk is fixed. In reality, the disk temperature cools 
significantly as M decreases, resulting in the inward migration 
of the various sublimation lines (see Garaud & Lin 2007, for 
instance). This will also affect the shape of the surface density 
profile in the inner disk; feedback between the disk temperature 
profile and the evolution of the surface density of grains will be 
the subject of future work. 

6. DISCUSSION 
6.1. Discussion of the particle size distribution function. 

The fundamental assumption underlying this work is that of 
the maintenance of a power-law particle size distribution func- 
tion at all radii, throughout the disk lifetime. The assumption is 
justified exactly only if the collisional timescale (note: not the 
growth timescale, which is naturally a factor of 1 /e larger) is 
shorter than the drift timescale for each size-bin. Whether this 
is in fact exactly true is certainly unlikely, but neither is it par- 
ticularly relevant. The correct questions that should be asked 
are: (i) how far from equation ([Sj is the true size-distribution 
function in a disk, and (ii) how do deviations from ([8]l impact 
the conclusions from this paper? 

Question (i) is a fundamental question, with implications 
reaching far beyond the scope of this paper. Attempts at an- 
swering it have come from various angles including both direct 
observations (in disks, but also in molecular clouds, in the ISM, 
as well as in our own solar system) and numerical experiments. 
As mentioned in 32.31 the observational evidence and theo- 
retical motivation for a power-law size-distribution function is 
strong but limited to more-or-less spatially isotropic and homo- 
geneous cases where there exist no systematic size-dependent 
drift or settling velocity which could act as a size-filter Numer- 
ical simulations of the coagulation-shattering balance in similar 
conditions also unanimously agree on the power-law distribu- 
tion, with indices varying slightly depending on the model as- 
sumptions but never straying too far from -3.5. 

Unfortunately, there has only been one study (Suttner & 
Yorke 2001) that self-consistently includes a complete param- 
eterization of the coagulation/shattering balance together with 
radial drift and vertical settling of the particles in an accretion 
dislsQ. Suttner & Yorke (2001) studied the formation of a pro- 
tostellar accretion disk through the collapse of a uniformly ro- 
tating molecular cloud core, and closely followed the evolu- 
tion of the grain size distribution function at every spatial po- 
sition throughout the collapse phase (first 10,000 years). They 
found that accretion shocks play an important role in limiting 
the growth of the grains; they also found, as expected, that 
the assumed sticking efficiency essentially governs the maxi- 
mum grain size achievable. The dust size distribution functions 
computed vary strongly with height above the disk: they show 
clear evidence that larger grains tend to be found in the mid- 
plane, while regions high above the mid-plane remain close to 
the MRN-derived initial conditions. This can be attributed to a 
combination of settling of the larger grains as well as preferen- 
tial in-situ growth. The mid-plane regions in their simulations 
appear to be largely depleted of small grains, which would be 
evidence for strong deviations from the power-law structure I 
assume. However, one may wonder whether this depletion is in- 
deed true in a real disk, since Suttner & Yorke neglect diffusion 

^ The study by DuUemond & Dominik (2005) does not include radial drift, and 



of the smallest grains by the gas turbulence (which could eas- 
ily bring small grains back towards the mid-plane from higher 
regions of the disk). 

In conclusion, one should bear in mind that the assump- 
tion made in selecting a power-law size distribution function 
is probably not always strictly justified (in particular for larger 
particles). But given the enormous computational advantage of 
this approach, it should be thought of as an acceptable trade-off 
between models in which the full coagulation/shattering equa- 
tion is solved, and models in which only one particle size (or a 
few particle sizes) are considered. 

Question (ii) can easily be answered by identifying where 
in the proposed model the assumption of a power-law size- 
distribution function is used. 

As mentioned in 32.51 the minimum particle size s^ia plays 
no role in the dynamical evolution of the solids in the disks, 
as long as ^ ^max, which is the case for the MRN size- 
distribution function, and therefore likely to continue being the 
case throughout the disk evolution. Thus whether the smallest 
grains are slightly depleted or not compared with the proposed 
power-law distribution function really does not matter. 

Collisional growth is essentially dominated by encounters 
between particles of similar sizes: even if collisions with 
smaller particles are more frequent, the mass gained is much 
smaller. This is easily seen mathematically in the derivation 
of the growth rates di^ax/df in 32.51 where the integral over 
all possible impactor sizes is always dominated by the largest 
particles, except possibly when 5fmax ^ 1, in which case the 
integral is dominated by particles of intermediate size. Another 
way of seeing this is that while the assumptions concerning the 
grain size distribution function used to derive equation |36] are 
very different from those of Stepinski & Valageas (1997) who 
assume that the distribution function is strongly peaked around 
a single-size, the outcome is the same whithin some factors of 
order unity. 

The systematic radial motion of decoupled particles (S'fmax ^ 
1) is the only place where the assumption made has a signifi- 
cant effect on the model results. If one assumes that all par- 
ticles have a single size i^ax, then Mp(imax) « l/5fn,ax while if 
one considers the mass-weighted average motion of all parti- 
cles following the proposed size-distribution function dHJ then 
Mp cx 1 / -\/5f,Tiax which can be significantly larger. As mentioned 
earlier, this accounts for the fact that when a continuum of par- 
ticle sizes is taken into account, intermediate-size grains (with 
St{s) ~ 1)) do rapidly drift leading to a non-neligible mass flux, 
even if the largest ones are fully decoupled. 

The scaling Mp cx 1 /V^fmax clearly depends on d«/ds cx s"^^; 
whether this exact scaling really applies to disks certainly is de- 
batable, but I would argue that the general picture of large bod- 
ies being eroded by collisions and leading to a non-negligible 
mass-flux even when the larger bodies themselves do not drift 
has to hold. The only caveat is that the collision rate becomes 
null when the bodies reached isolation mass; thus one should, 
for self-consistency, set Mp — + in this limit, which was not done 
here (otherwise, smaller particles coming from larger radii arti- 
ficially accumulate in the inner regions). A way forward would 
be to combine the model proposed here with an N-body code, 
in which particles are treated using a size-distribution until they 
reach embryo size, then taken out of the distribution and indi- 
vidually followed using N-body simulations. This could be the 
subject of future work. 

treats shattering in a very simplistic way. 
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6.2. Heavy element retention efficiency in the UV-switch 

model 

The simple equation for the evolution of the total mass of 
solids JTSI ) can be used to derive the final contents of the disk 
after complete gas dispersal when caused by photo-evaporation 
from the central star. 

If grains in the outer disk remain fully coupled to the gas 
throughout the lifetime of the disk (i.e. if tp ^ Tgap, or equiv- 
alently, ^maxO <C 1/im) then the amount of material left after 
complete dispersal of the gas can easily be estimated by the 
amount of solids left at f = Tgap (see AA07). Therefore an or- 
der of magnitude estimate for the heavy-element retention effi- 
ciency is simply 

Mp(Tgap) ^/ Mo ^^g^ 

Mpo V2tvM„/ 

where Mpo = ZqMq. For a fixed initial disk mass, this estimate 
depends weakly on Ro; this could explain the very low dis- 
persion observed in the heavy element retention efficiency of 
evolved systems (Wilden et al. 2002) despite the vastly differ- 
ent disks dispersal timescale required by SED observations. 

To refine this estimate and quantify the effect of the ini- 
tial grain size distribution on the heavy-element retention ef- 
ficiency, I integrate ( |75] l from f = to f = Tgap for a wide variety 
of initial conditions (Mq, Rq). The results are shown in Figure 
[T2I Four cases are considered, with varying initial parti- 
cle size .Smaxo and turbulent at. The weak dependence on the 
initial conditions of the disk (Mq, Rq) for given values of at 
and SmaxO suggested by equation ( fTSl l is naturally still present: 
it appears that even with Mo and /?o varying by two orders of 
magnitude, the remaining amount of solids does not vary more 
a factor of a few. In the fiducial model for instance (5,^3x0 = 1 /^m 
and at = 0.01) the typical amount of solids left is of the order 
of a few Earth masses for most plausible values of Mq and R(). 

As suggested by the results of ^5.3.21 the retention efficiency 
drops dramatically for larger initial grain sizes, but naturally 
saturates near the value given by equation ( fTST l for very small 
grain sizes (e.g. compare the results for imaxo = 3/im with the 
results for imaxO = 0.3/im). Also suggested by ( fTSl l is the de- 
pendence of the phenomena on at: for the smaller value of 
at = 0.001, one could expect up to a few tens of Earth masses 
to be left behind in the disk, a value that begins to be consis- 
tent with the amount of solids left in the Minimum Solar Neb- 
ula model augmented with the mass of the Oort cloud. Thus it 
appears that consistency of this idea with our own solar system 
would strongly favour a model with imaxO is no larger than 1 /im, 
and at is preferably of the order of 0.001. 

How reliable is this estimate given the simplifications made 
in the model? Conveniently, given that the total disk mass re- 
sides mostly in the outer disk, and that the amount of mate- 
rial remaining in the inner disk is ultimately controlled by the 
mass flux from the outer reservoir, the mass estimates given 
here are reasonably independent of the physics of the inner 
disk (including sublimation/condensation, but also dead zones, 
etc.). There is however an important caveat; as mentioned ear- 
lier, large protoplanetary embryos which have reached isolation 
mass are fully decoupled from the disk dynamics (both in terms 
of their drift velocity and in terms of their collision frequency), 
and are not well-modeled by the size distribution function pro- 
posed. The mass contained in these embryos could add a few 
Earth masses to what is presently estimated, but only in the 
inner disk. The total mass of solids which ends up being shep- 



herded out to the outer solar system is not affected by this prob- 
lem, and is therefore reliably estimated by this method. 

7. CONCLUSIONS 

This paper presents a new algorithm modeling the evolution 
of gas and solids in protostellar, as well as some reliable quan- 
titative analytical estimates for the outcome of the numerical 
simulations. 

The global disk evolution paradigm is well-reproduced by 
the numerical solutions. Well-known results are recovered, 
such as the two-timescale gas evolution (Clarke, Sotomayor & 
Gendrin 2001; Alexander, Clarke & Pringle 2006a), the rapid 
growth of solids in the inner disk (Suttner & Yorke, 2001, 
Stepinski & Valageas, 1997, DuUemond & Dominik, 2005), 
the shepherding of smaller particles by the retreating hole front 
(Alexander & Armitage 2007), and the accumulation of ma- 
terial near the sublimation lines (Stevenson & Lunine, 1988, 
Stepinski & Valageas 1997, Ciesla & Cuzzi, 2006). 

Novel conclusions of this paper are: 

(i) The evolution of the mass of solids in the disk is essen- 
tially controlled by a reservoir of small grains at large radii. A 
well-tested empirical formula for the total solid mass in the disk 
is given in equation dTSl l. 

(ii) The heavy-element retention efficiency after gas dispersal 
is controlled by the remaining amount of solids left at the time 
of gap opening, and is found to vary weakly with initial disk 
conditions but very sensitively with initial particle size ^maxO- 
The remaining amount of solids in the Minimum Solar Neb- 
ula combined with the mass of the Oort cloud is inconsistent 
with iiTiaxO greater than 1/im, and would tend to prefer a value 
of at ~ 0.001. 

(iii) The strong dependence of the gas dispersal timescale on 
the initial mass and radius of the disk combined with the weak 
dependence of the heavy element retention efficiency on the 
same parameters could simultaneously explain the wide diver- 
sity of the SED observations with the very low dispersion of 
the stellar metallicities observed in the Pleiades (Wilden et al. 
2002). 

(iv) Rapid grain growth in the inner disk implies that solids 
in the form of large planetesimals (with sizes ranging from a 
few meters to 1000 km) are locally retained after gas dispersal. 
The presence of a population of large planetesimals in the inner 
regions of debris disks has been inferred from interferometric 
observations of the presence of dust despite its very short radi- 
ation blowout time (specifically in TW Hya, by Eisner, Chiang 
& Hillenbrand 2006) 

(v) All the small grains are swept by the retreating gas front 
at the edge of the hole and shepeherded outward. When the ac- 
cumulated surface density of grains approaches that of the gas, 
the gas becomes unable to continue moving the grains and will 
most likely leave them behind. This could explain the system- 
atic formation of narrow dust rings at large radii. 

(vi) Regulation of the solid mass flux by the "leaky reservoir" 
implies that any local surface density enhancement (or "peaks") 
near sublimation lines must be accumulated very early on in 
the disk lifetime, more precisely during the initial phase when 
the disk dynamics are still out of equilibrium (~ first 10,000- 
100,000 yr). The gradual decay of the inner disk gas density 
through photo-evaporation could then trigger gravitational in- 
stabilities at well-separated times as each peak respectively ap- 
proaches unit metallicity. 

(vii) Possible constraints on the sticking efficiency of parti- 



21 




1 10 100 1 10 100 

Ro (AU) Ro (AU) 



Fig. 12. — Heavy element retention of disks after total photo-evaporation of the gas as a function of the initial conditions of the disk, for various initial particle 
sizes and turbulent parameter at. The color scheme is the same for aU four plots. The soUd lines (at 1, 2, 5, 10, 20 and 50%) mark the retention efficiency, namely the 
percentage of heavy elements remaining in the disk compared with its initial content ZqMo. The dotted lines (at 10"*, 2 x 10"*, 5 x 10"*, 10"', 2 x 10"', 5 x 10"', 
10"^ and finally 2 x 10"*M«) foUow the colored contours and mark the actual total mass of the remaining solids. Note that IM© ~ 3 x 10"*Mq so that the 10"^ 
contour corresponds to about 33 M^. 



cles in the turbulent conditions found in protostellar disks could 
be derived from spatially resolved observations of the grain sur- 
face density of nearby disks (see Appendix C). 

Direct comparison of the model predictions with disk SEDs 
is the subject of Paper II. 
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APPENDIX 

APPENDIX A: PHOTO-EVAPORATION MODEL 

Alexander & Armitage (2007) studied the following model for the mass loss rate from photo-evaporation from the central star, 
based on the works of HoUenbach et al. (1994), Font et al. (2004), as well as Alexander, Clarke & Pringle (2006a and 2006b). 

Prior to the removal of the gas in the irmer disk, the mass loss rate is caused by the diffuse UV field reprocessed high in the disk 
atmosphere. It is equal to 

t^{r) = 2noir)uiir)i^mH , (1) 
where /i is the mean molecular weight of the ionized gas (taken to be /i = 1.35), oth = 1.67 x 10~'^*g is the mass of the Hydrogen 
atom, the column density at the base of the ionized layer no is taken to be (see Font et al. 2004) 



(r//?g)15/2 + (r//;g)25/2 



no(»-) = 0.14| 



47raB/?2 / 
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with $ is the photo-ionizing flux (in photons per second), = 2.6 x 10~'-'cm^/s, Rg = GM^:/cj and c, is the sound speed of the 
ionized gas (taken to be lOkm/s in the fiducial model). The launch velocity ui(r) is taken to be (see Font et al. 2004) 

0.2457 



Mi(r) = 0.3423c,exp 
= Qifr<Q.lR 



-0.3612 



if r >0.1«„ 



. (3) 
After clearing of the inner disk, the mass loss rate is mostly caused by direct photo-ionization of the inner edge of the remaining 
disk, and accordingly changes to 

-1-1 / * ^ 1/2 //,\-l/2 / ^ \-2-42 

(4) 



Sw('',f) = 0.47 /imnc/ 



1 +exp 



-'-hole(f) 



47raBri;„i^(f) 



'"hole(0 



where rhoie(f) is the hole radius, well approximated to be the radius for which S(r,f) drops below 10 ^g/cm^ (Alexander, private 
communication). 



APPENDIX b: numerical method 

Equations ([T]i, © (for each species), ( [36] l or ( l44b . and finally ( |46] l (for each species) are evolved simultaneously in time using the 
following approach. 

The independent variables r and t are first normalized to 1 AU and to Tau = 27r/fiK( 1 AU) respectively. Following Pringle, Verbunt 
& Wade (1986), a new variable is introduced 

1/2 /IN 

y = '•au ' (1) 
upon which space is uniformly discretized. Next, new dependent variables G, G[, and GJ, are defined as 



G = 5] and g; =/«-^I] 



v,p 



With these modifications, equations ([Hi and ( l46b can be rewritten as 

dG 



dt 



AU 



(2) 



(3) 



(and similarly for each vapor-form species) where t^u = f /7au- The evolution equation for the solids becomes 



dtAv dy 



TT -2 _i 91nG iTTl^J 

■^a,^hAv(Sc,fi-3I)— 



1 9G' 



Gi 



(4) 



for each solid-form species. This formulation emphasizes the pure diffusive terms, which cause numerical instabilities if not treated 
adequately. For simplicity, / = li^/lnSt^) and similarly for 7. 

Equations ( |36l ) or (|44| depending on the regime considered involve no spatial derivative, and are therefore simply integrated 
forward in time using a simple explicit second-order Adams-Bashforth scheme. Equations dUi and (O are integrated forward in 
time using a centered implicit scheme for the diffusion terms, an upwind first-order scheme for the advection terms and a second- 
order Adams-Bashforth scheme for all other terms. With this particular semi-implicit scheme, it is found that using about 1000-4000 
meshpoints (depending on the radial resolution required) and a typical timesteps Af ~ 0. 1 - lOyr yields a stable and accurate solution. 

APPENDIX C: A QUICK-AND-DIRTY ESTIMATE FOR Sp(r,f) 

To get an estimate for T,p{r,t), one could use equation ( ISOl l together with equation ( |76] |. Unfortunately, the value of 5fmax('",0 
is unknown since it depends both on J'max('",f) and on S(r,r), neither of which are particularly well approximated in inner and 
intermediate regions of the disk. Nonetheless, let's construct a low-quality estimate of 5'fmax('", using the following algorithm: at 
each radius r and time f , 

2 



Ant) 



mi,so(r,f) 



max 



lyear 



AU 



Ml. 

Mq 



27rr^ZoS(r,f)fe 



2 
3M^ 



1/3' 



-1 3/2 



SmAr,t) := min 



3OTiso(r,f) 
47r/9s 



1/3- 



(1) 



/27r7E(r,f) 

where miso{r,t) is an estimate of the local isolation mass, and where S(r, f) is given by the self-similar solution ( l55l ). The resulting 
predicted solution for Sp(r,/) is shown in Figure[T3]at t = 4.5Myr, for the three initial particles sizes considered in the text (1,3 and 
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Fig. 13. — Trae numerical solution (solid lines) and approximate estimates for the total suiface density of solids, using the algorithm shown in Appendix C. The 
estimate for l^m size particles is shown as a dotted line, 3/^m as dot-dash line and lOfim as a dashed line. Note that some of the global trends (in particular the 
overall normalization of the curve, and the surface density dip) are well-reproduced by the estimate but that quantitative agreement is poor 

10/zm). While the global features of the solution are fairly well reproduced, the quantitative predictions are clearly off by up to an 
order of magnitude in regions I and II of the disk, and even larger in region III where equation ( l76b looses validity. Note that one of 
the features that is well-reproduced by the solution is the sharp decrease in the surface density of particles in the intermediate regions 
of the disk. This feature corresponds to radii where 5ftnax — 1, which is approximately where 



0.01 



0.01 



0.01 



IMyr 



Interestingly, raip is entirely independent on the initial conditions of the disk (Mq and Rq). Since Zq, M-^ and t should be fairly well 
observationally determined by the stellar properties, possible measurements of rdip may help constrain the product e^at (at least 
within an order of magnitude). 
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